Axially symmetric focusing as a cuspoid diffraction catastrophe: 
Scalar and vector cases and comparison with the theory of Mie 
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An analytical description of arbitrary strongly aberrated axially symmetric focusing is developed. 
This is done by matching the solution of geometrical optics with a wave pattern which is universal for 
the underlying ray structure. The corresponding canonical integral is the Bessoid integral, which is a 
three-dimensional generalization of the Pearcey integral that approximates the field near an arbitrary 
two-dimensional cusp. We first develop the description for scalar fields and then generalize it to the 
vector case. As a practical example the formalism is applied to the focusing of light by transparent 
dielectric spheres with a few wavelengths in diameter. The results demonstrate good agreement with 
the Mie theory down to Mie parameters of about 30. Compact analytical expressions are derived 
for the intensity on the axis and the position of the diffraction focus both for the general case and 
for the focusing by microspheres. The high intensity region is narrower than for an ideal lens of 
the same aperture at the expense of longitudinal localization and has a polarization dependent fine 
structure, which can be explained quantitatively. The results are relevant for aerosol and colloid 
science where natural light focusing occurs and can be used in laser micro- and nano-processing of 
materials. 

PACS numbers: 42.15.Dp, 41.20.Jb, 42.25.Fx, 81.16.-c 



I. INTRODUCTION 

Axially symmetric focusing of wave fields occurs in var- 
ious areas of science, since physical systems often possess 
an intrinsic rotational symmetry. In particular, the elec- 
tromagnetic field enhancement by small spherical parti- 
cles is important in many situations. Spheres have mini- 
mal surface energy for a given volume and thus are natu- 
rally formed as a result of phase separation, for example 
as aerosols or colloids. Applications of colloidal micro- 
spheres in photonic crystals and photonic crystal slabs 
led to an explosion of the experimental and theoreti- 
cal studies of their optical properties. ^'^ The majority of 
these investigations concentrate on their collective prop- 
erties in a periodic arrangement. Single microspheres 
are used as high quality optical resonators and as agents 
that allow controlled and highly localized wavelength- 
dependent field enhancement for non-linear optical stud- 
ies and in resonance spectroscopy^!^ Here, the emphasis 
is placed on the cigenmode analysis and the distribution 
of the field within the sphere or in the immediate vicinity 
of its surface. 

Lately, it was demonstrated that self-assembling ar- 
rays of transparent colloidal microspheres can be em- 
ployed for high-throughput laser-assisted micro- and 
nano-structuring of materials.^'^'^ Similar effects were ob- 
served in experiments on dry laser cleaning, where such 
particles are used as controlled contaminantsiSiSiiS This 
necessitates better understanding of the focusing of light 
by microspheres with diameters of several wavelengths. 
Only few rigorous results are available for the interme- 
diate range of sphere sizes and distances from the par- 



ticle. The majority of analytical approaches either deal 
with the properties of the eigenmodes, or refer to the in- 
tegral characteristics and/or to the far field behavior lii 
Mie resonances were analyzed on the basis of advanced 
geometrical optics^^ and detailed numerical calculations 
for transparent spheres of several wavelengths in size were 
performed in connection with the use of laser tweezers in 
biology^^ and for needs of aerosol sciencesi^ 

In this work we develop a theoretical description for 
an arbitrary non-paraxial strongly aberrated axially sym- 
metric focusing and apply it to the case of dielectric mi- 
crospheres. Our emphasis is on the fine structure of the 
field distribution in the exterior of the sphere up to the 
focal region, which can be used to control and improve 
the concentration of energy. 

Strong spherical aberration makes the focusing non- 
trivial. Usually, the exact solution is obtained using the 
Mie theory, which does not give much of a physical 
insight as it requires the summation of a large number of 
terms in a multipole expansion even for moderate sphere 
sizes. At the same time, the main focusing properties of 
transparent dielectric microspheres originate rather from 
the picture of geometrical optics. 

One might think that in the lowest approximation a 
small sphere acts as an ideal lens. However, in the range 
of sizes we are interested in, this picture does not even 
provide a description which is qualitatively correct. Also 
classical formulas for weak spherical aberration^^ do not 
yield useful results for the field behind a sphere: They 
predict that the maximum intensity is kept unchanged 
and its position does not depend on the wavelength. 

Our approach, following the method of uniform cans- 
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tic asymptoticsf"- is based on the canonical integral for 
the cuspoid ray topology of strong spherical aberration. 
Though this Bessoid integral — a member of the hier- 



archy of diffraction catastrophe 



.18.19 



appears natu- 



rally in the paraxial approximation, it can be used to de- 
scribe arbitrary axially symmetric strong spherical aber- 
ration by appropriate coordinate and amplitude transfor- 
mations. For angularly dependent vectorial amplitudes 
the formalism uses higher-order Bessoid integrals. 

The Bessoid integral is the axially symmetric general- 
ization of the Pearcey integral,^" which plays an impor- 
tant role in many short wavelength phenomenal^ There- 
fore, the present approach can be applied in various areas 
of physics where axially symmetric focusing is of impor- 
tance, e.g., acoustics, semiclassical quantum mechanics, 
radio wave propagation and scattering theory. 



II. THE BESSOID INTEGRAL 

A. Definition 

We first consider the diffraction of a scalar spherically 
aberrated wave on a circular aperture with radius a in 
the plane z = —f around the z-axis, where / is the focal 
distance. The origin of the coordinate system is put into 
the focus F. In cylindrical coordinates (p, z), the parax- 
ial Fresnel-Kirchhoff diffraction integral^^ yields the field 
amplitude 



ikUo 



J \ ~)^ Pid/Oi. (1) 



Here Uq is the amplitude of the incident wave in the cen- 
ter of the aperture, k is the wavenumber (fc — 27r/A, 
where A is the wavelength) and pi is the distance from 
the axis on the aperture. The Bessel function Jq comes 
from the integration over the polar angle ip. The param- 
eter B in the exponent determines the strength of the 
spherical aberration. For B > the diffraction focus 
shifts towards the aperture, while B = corresponds to 
ideal focusingii^ 

We i ntroduce the dimensionless coordinates pi = 
V4JBpi,R = ^/WJIB p/f and Z = y^k/lBz/f and 
consider an infinitely large aperture. Then the field 
becomes proportional to the Bessoid integral^ 
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FIG. 1: Absolute square of the Pearcey integral Ip (top) and 
the Bessoid integral / (bottom). The latter is proportional to 
the field of a spherically aberrated wave within small angles 
approximation. 



Its absolute square is shown in figure ^ In the Carte- 
sian representation xi — pi cos (p and yi — pi sin (p are 
dimensionless coordinates in the plane of integration. Ex- 
pression ||2J) is the axially symmetric generalization of the 
Pearcey integral^S 



/p(X, Z) 



1 



/2 7r 



Xxx + Z ^ + ^ 



dxi, (5) 



which is also shown in figure ^ 

Both integrals correspond to so-called diffraction 
catastrophesM^i^ii^ Their field distribution contains 
caustic zones where the intensity predicted by geomet- 
rical optics goes to infinity. The Pearcey integral corre- 
sponds to a cusp caustic, i.e., a single one-dimensional 
curve in a two-dimensional space, and does not reveal 
a high intensity along the axis, while the Bessoid inte- 
gral corresponds to a cuspoid caustic, i.e., to a surface 
of revolution of the cusp in three dimensions, as well as 
the caustic line up to the focus F eA z = Z = Q. The 
equation of the cusp is given by the semicubic parabola 



27 + 4 Z'^ = . 



(6) 



Henceforth we will apply the term cusp also for the whole 
cuspoid. A caustic is denoted as stable, if it does not 
change its topology under small perturbations. This is 
the case for the Pearcey integral. The Bessoid integral 
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FIG. 2: (a) 3-ray region inside the cuspoid (dashed line), (b) 
1-ray region outside. The z-axis is represented by a dashed- 
dotted line. 



corresponds to a structurally unstable caustic, because 
an infinitely small perturbation will destroy the radial 
symmetry and the axis will not be a caustic zone any 
longer. It is, however, stable on the class of axially sym- 
metric wavefronts. 

The cusp is the envelope of the family of rays. The 
latter correspond to the points of stationary phase in the 
Bessoid integral, i.e., those points where the two first 
partial derivatives with respect to R and Z of the phase 
(/) in © vanish. Inside the cusp, for 27 R"^ + AZ^ < 0, 
three rays (tangents to the cusp) arrive at each point of 
observation P = {p,z), and outside, for 27 R'^ + AZ^ > 
0, there is only one real ray (figure Thus, the cusp 
forms the border between the lit region and the (partial) 
geometrical shadow, where two rays merge. 

Without loss of generality, we assume that all rays lie 
in the meridional plane ip ^ (yi = 0) and hence corre- 
spond to the roots xij {j — 1, 2, 3) of the cubic equation 



stationary laa 

™ n i -/"j + i f sign Hj 

where the summation runs over all real rays, i.e., m = 1 
(lit region) or m = 3 (shadow). The phase (j)j is obtained 
by inserting the j-th stationary point {xij,yij — 0) into 
the l@J. The determinant and signature of the Hessian 
are given by 

det Hj = + 4 xl j Z + 3 xfj , (9) 
signHj = sgn(-Z - 3x1 j) + sgn(-Z - xl ^) . (10) 

Near the cusp the Bessoid integral shows an Airy-type 
behavior typical for caustics where two real rays disap- 
pear and become complex. 

One can also derive a different approximation valid on 
and near the caustic axis, i.e., for Z < and small R 
(appendix : 

I{R,Z)^^MRV^)e^^ erfc(^-e'4j. (11) 

Here erfc is the complementary error function, which 
can also be written in terms of Fresnel sine and cosine 
functionsj22i Expression (|11|) becomes exact at the axis 
R = 0, where Jo(0) = 1. It shows that near the axis 
the Bessoid integral is virtually a Bessel beam^S, with a 
variable cross section. 



C. Numerical evaluation 



R + Zxi+xl = 0, (7) 

which are given by Cardan's formulas^ On the axis, 
i? = 0, a cone formed by an infinite number of rays con- 
verges. These rays originate from the circle x^+yf = —Z 
on the aperture. They are all in phase and produce a 
high intensity along the axis (compare the two pictures 
in figure ^ . The oscillations occur due to interference 
with the ray propagating along the z-axis. One can also 
directly observe that the Bessoid integral has the topol- 
ogy of spherical aberration as the maximum of intensity 
does not lie in the geometrical focus Z = but is spheri- 
cally aberrated to a negative value of Z, i.e., towards the 
aperture. 



B. Asymptotic expressions 

Off the caustic — away from the cusp and the focal 
line — the Bessoid integral (PJ can be approximated by 
the method of stationary phase. As the integrand in Q 
is highly oscillatory, the only significant contributions to 
the integral come from those regions where the phase is 



As the Bessoid integrand is highly oscillatory, its eval- 
uation for the whole range of coordinates R and Z is 
non-trivial and of large practical importance. Direct nu- 
merical integration along the real axis and the method 
of steepest descent in the complex plane both have their 
disadvantages. By far the fastest technique is based on 
the numerical solution of the ordinary differential equa- 
tion (derivation in appendix IB|! ^^'^^ 

Lr- Z Ir + iRI = 0. (12) 

Indices denote (partial) derivatives and L = Irr + Ir/R 
is an abbreviation for the radial Laplacian applied onto 
/. The three initial conditions at i? = are 

/(0,Z) = ^e'^erfc(^-e'4j, (13) 

/i?(0,Z) = 0, (14) 
L(0,Z) = Z/(0,Z) + i. (15) 

I{0,Z) was taken from ifTT)) . Ir{0,Z) vanishes due to 
symmetry, and the last condition arises from the fact 
that the Bessoid integral satisfies the paraxial Helmholtz 
equation 2ilz + L = 0, where Iz is calculated from l|13|l . 
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In the literature the Pearcey integral was calcu- 
lated by solving differential equations^S by a series 
representationSli and by the first terms of its asymp- 
totic expansionii^ The Bessoid integral was expressed in 
terms of parabolic cylinder functions''^ and as a seriesiS^ 
The latter work gives reference to an unpublished work 
of Pearcey^Si stating that differential equations for the 
Bessoid integral were employed there. 



D. Geometrical optics for the cuspoid 

In geometrical optics, the rays carry the information 
of amplitude and phase. The total field in a point P is 
given by the sum of all ray fields there. A ray's field at 
P is determined by'^^ 



U{P) = Uo 



7T 



(16) 



where Uq is the amplitude at some initial wavefront, tp 
is the eikonal, and J is the generalized geometrical di- 
vergence, which can be calculated from flux conservation 
along the ray. For a homogeneous medium with constant 
refractive indes 
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J = 



(17) 



Rm, Rs are the main radii of curvature at the point P 
and i?mOi Rso are the radii on the initial wavefront, where 
U = Uo. 

When a ray touches a caustic, its radius of curvature 
(the geometrical divergence in the general case) changes 
the sign and the ray undergoes a caustic phase delay^^^^ 
of — 7r/2, which is taken into account by the proper choice 
of the square root in ((T?^ . When a ray touches several 
caustics, these delays must be added. The total caustic 
phase shift, denoted as A(p, can be explicitly written in 
the phase. For the cuspoid topology and ray numbering 
(j = 1,2,3) according to figure 13 we obtain: 



with 



U{P) = C/o ■ 



^ i /c -0 

77 



(7o 



^ i i/f + i Ai^ 



— vr for J = 1, 
for j = 2, 
-7r/2 for j = 3. 



(18) 



(19) 



Ray 1 touched the cuspoid and the focal line, ray 2 is not 
shifted, and ray 3 touched the cuspoid. 



III. RELATION BETWEEN GEOMETRICAL 
AND WAVE OPTICS 

A. Matching with the Bessoid integral 

If we have found the phases if = kip and divergences 
J of the rays, the (scalar) geometrical optics solution 



with an axially symmetric 3-ray cuspoid topology can be 
written as 



U{v) 



E 



Uq^j e 



(20) 



Here r = (p, z) are the real-space coordinates and we 
have allowed for different initial amplitudes Uqj of the 
rays. This field shows singularities at the caustic, espe- 
cially on the axis, which is the most interesting region for 
applications. 

We want to describe arbitrary axially symmetric fo- 
cusing by matching the solution of geometrical optics 
(where it is correct) with a wave field constructed from 
the Bessoid integral ©, which naturally appears in the 
paraxial approximation and is finite everywhere, and its 
partial derivatives and Iz {method of uniform caustic 
asymptotics) . We make the Ansatz^ 



U 



AI 



1 



■AnL 



1 



R Jfl 



Aziz 



(21) 



The yet unknown arguments of the Bessoid integral and 
its derivatives are R = (i?(r), Z(r)). A{y), Afl(r) and 
Az{v) are three amplitude factors and x(r) is a phase 
function. (The indices R and Z in the amplitudes do 
not indicate derivatives.) Now the geometrical optics 
solution (|20|l is matched with the stationary phase ap- 
proximation of l|21|l by equating the amplitudes and 
phasesilMi 



A{v) + AR{r) 0k(R, t,) + Az{r) 0z(R, tj 



V'j-W =x(r)+0(R,t,). 
bR and 4>z are the partial derivatives of I0J and 

I (,14 sign Hj 



(22) 
(23) 



(24) 



where the determinant and signature of the Hessian are 
written in and (|10(l . respectively. Outside the cusp, 
the rays 2 and 3 are complex and the general definition 
of Hj is more subtle, namely 



020, i Y 4^02,3 



(25) 



with = d'^(f)/dx\, 002 = d^^jdyX (the index j de- 
notes substitution of the j-th point of stationary phase 
as argument). 

The three points of stationary phase were denoted as 
tj = (tj, 0), where the tj are given by the (correctly or- 
dered) Cardan's solutions of iQ, i.e., of 

i? + Z< + f''=0. (26) 
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Note that they are functions of the Bessoid coordinates, 
tj — tj (R) , and the latter depend on the real space coor- 
dinates: R = R(r). The partial derivatives with respect 
to R and Z in H22I) must be evaluated in such a way as 
the tj were held constant, although they are functions 
of R themselves. The conditions H22|) and H2^^|l give 6 
equations for the 6 unknowns R, Z, A, Aji, and Az- 
It is convenient to solve (|23|l . that is 



fj =x- 



Rtj--Zt] 



(27) 



using quantities that are permutationally invariant with 
respect to the roots tj i ^^i?^ This yields 



""-'^-A-JZ 



Z^±\l\ \/-2 sgn(63) ^b^q + 2Vd, 



(28) 



D 



262-9 
1 



^2 9 + 



6 



z^ 



where sgn(Z) = sgn(Z4 - 2462)- The bi (/ = 1,2,3) are 
given by 61 = il/3)Y.^.^^ipj , 62 = " bi)^ and 

63 = J2'j=ii'Pj ~ bi)^- The quantity q (sometimes called 
discriminant) can be expressed in different ways: 



6 b. 



3^3 



— i?^ {27R^ + 4:Z-^) 



-2 [ipi - LP2T {^2 - ipaf (<^3 - ViY 



(29) 



Hence, it vanishes exactly at the caustic where two phases 
are equal. At the cuspoid ip2 = ^3 (27 R"^ + iZ^ = 0) 
and on the axis ipi ~ ^3 {R — 0). 
The solutions of if^ . that is 



A - tj Ar 



\t]Az 



(30) 



ari 



,35 



Ar^ 



A, 




t2 ts 



Jl {t3-ti){ti 

h; ti 



t2) 



2Uo 



Jl (i3-ti)(ti 
Th 1 



t2) 



(31) 



Jl it3 - h) [tl - t2) 



where the cyclic terms permutate the numbering of rays: 
(1, 2, 3) -> (2, 3, 1) (3, 1, 2). The Bessoid matching so- 
lution H21|l does not show the divergences of geometrical 
optics. 

Note that this method utilizes also the so-called com- 
plex rays which have less apparent physical meaning. It 
turns out that both real and complex rays provide the 
geometrical skeleton for the wave fleshi^ 



B. Expressions on and near the axis in the general 

case 

All formulas can be strongly simplified on and near 
the axis inside the cuspoid (small p, z < 0) . The Bessoid 
coordinates have the simple form (appendix Q 



R 



k p sin P 



V2^iipi+ip3)/2 - ^2 
Z w -2 



•fi + ^3 



IP2 ~ -2 ^(^„p - <y9p , 



(32) 
(33) 



where /3 > is the local angle of the non-paraxial cone 
of rays with the axis and iy9np and denote the phases 
of the non-paraxial rays and the (par) axial ray, respec- 
tively (see figure El in appendix . The simple natural 
combination 



k p sin (3 



(34) 



also appears in the near axis approximation for the 
Bessoid integral (|llll . On the axis [p — Lpi — (^3) 
we obtain i? = and Z — —2 y/^pi — (^2- 

The results (|32|l - (|34|) have transparent physical mean- 
ing. Indeed, near the axis the largest contribution to 
the field comes from the converging cone of non-paraxial 
rays (similar to ray 1) that intersect the axis at an angle 
/3. If the angle (3 is constant and all rays have the same 
intensity, the result is the Bessel beam.'^^ Such beams 
have a propagation constant along the z-direction equal 
to the z-component of the wavevector of the plane waves 
which form them and correspondingly the argument of 
the Bessel function (cylindrical analog of a plane wave) 
is equal to k p sin (3. As the angle (3 gradually changes 
for the spherically aberrated wave, so does the argument 
of the Bessel function. 

Additionally, there exists the axial ray 2, which is not 
present in the canonical Bessel beam (though it often 
appears in real experimental situations). The interfer- 
ence of this beam with the converging ray cone results 
in the intensity oscillations along the axis (figure ^ bot- 
tom). Clearly, these oscillations are largely due to the 
phase difference ipnp ~ ^p- At large negative Z in 
erfc[| exp(if )] -> 2 [l-^^i^- exp(i^^=^)], and the oscil- 
lating behavior is governed by the phase of the exponent, 
which is equal to 3 7r/4 — ((/Jnp — <Pp). This clarifies the 
origin of expression H33|l . as it is Z"^ entering the final 
formulas. 

In particular, the global maximum is expected on axis 
at the first constructive interference of the axial and the 
non-paraxial rays. Because Z < in this region, the two 
terms of the erfc expansion are first in phase when the 
phase difference is fi~^p2 — 3 7r/4. The geometrical 
meaning of this result is that rays 1 and 3 are shifted by 
— 7r/2 as they touch the cusp. In addition, they acquire a 
further shift of — 7r/2 when crossing the focal line. But ex- 
actly on the axis only half of this delay has occurred yet. 
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which yields the 3 7r/4 difference. The numerical maxi- 
mum of the Bessoid intensity (absolute square) occurs at 
~ —3.051 and hence this yields the condition 



^ w 2.327, 



(35) 



which is close to 3 7r/4 « 2.356. 

The width of the focal line caustic, p^, is defined by 
the first zero wq ~ 2.405 of the Bessel function in 
Hence, with (jSH), 



Wn A 

° - 0.383 



k sin/3 



sin/3 



(36) 



In the geometrical optics picture the first minimum oc- 
curs when rays 1 and 3 interfere destructively, i.e., when 
their phase difference becomes it. This results in (pi — 
ip^ = it + it/2^ where the term tt/2 takes into account the 
caustic phase shift of ray 1: ^ {fi — 'Pa) /2 k sin f3 = 
0.375 A/ sin/3. Note, that this is smaller than the Airy 
spot for the same aperture anglei^ and large angles /3 
are indeed realized, e.g., in the case of the sphere studied 
below. 

Finally, we present an expression for the field 1211) on 
the axis. The equations for the amplitudes H31(l simplify 
tremendously f appendix ID|) and result in 



U = 



Uq^i V2 fcp sin/3 



i/-l 
Z 



0,2 



/J2 



(37) 



The structure of expression (|37|l helps to understand its 
physical meaning. It details the contribution of the cone 
of non-par axial rays, represented by ray 1 (first term), 
and the axial ray 2 (second term) to the overall structure 
of the field. Note that in the general case not only the 
angle /3, but also the amplitude of the converging cone 
may vary along z (Z), thus slowly modifying the proper- 
ties of the Bessel beam in the axial region. This enters 
(|37|l via amplitude transformations and is manifested by 
the presence of the initial ray amplitudes in both terms. 
Inside the cusp on the axis {z, Z < and p,R ^ 0) 
both \l\f7i and the ratio yfpl\fJi remain finite as the 
divergence of the paraxial ray 2 is non-singular, while the 
sagittal divergence Ji of the cone of non-paraxial rays 1 
is proportional to p. Due to the Bessoid matching pro- 
cedure the singularity of the converging cone is removed 
by the compensating factor y^. Along the axis the last 
term in (|37|l partly cancels with the second term in the 
parentheses of the first term. As a result, the on axis field 
behavior up to the focus is dominated by a single term 
proportional to the Bessoid integral /, which justifies the 
maximum condition (I35II discussed above. 



C. Angular dependences and vectorial problems: 
Higher-order Bessoid matching 

Often — especially in vectorial problems — there ex- 
ists axial symmetry with respect to the wavefronts, ray 



phases and generalized divergences, but not with respect 
to the amplitudes. In this case, new functions are re- 
quired to represent arbitrary angular dependence of the 
field. The natural generalization of (O are the higher- 
order Bessoid integrals^ with the non-negative integer 



JQ 



dpi 



(38) 

where Iq = I and J,„ are higher-order Bessel functions. 
The higher-order Bessoid integrals obey the recurrence 
relation 



T J I ™ 

— —J^m.R + — 
K 



(39) 



The integral is canonical for angular depen- 

dent geometrical field components U'^"^\p, z) slumps or 
U^-^\p,z) cos mip. In matching similar to H21() . 



(40) 

the angular dependence cancels. Here Am, A^r and 
AmZ are the higher-order amplitude factors, whereas 
Im,Rrr, and Im.Zm are partial derivatives of the higher- 
order Bessoid integrals Im- Since the latter can be writ- 
ten in terms of /qi it can be shown that the points of 
stationary phase, the matching of phases and thus the 
higher-order coordinates [Rm, Zm) and phases (xm) are 
identical with the original ones: 



Rm — R , Zm ^ Z , Xr, 



(41) 



From the physical point of view, this reflects the con- 
servation of the wavefront and thus the ray phases and 
divergences. 

The equations for the amplitudes have to be general- 
ized. The higher-order amplitudes Am, AmR and AmZ 
have the same form as (|31|l , but with an additional factor 
(iij)™ in each denominator, i.e., 



Am — — Uf)^i 
AmR — Uo^i — 
AmZ — 2 C/o,l ■ 



t2 ts 



(i<i)™(i3-ii)(ii-<2) 

t_i 

(iti)"(t3-ii)(ti -^2) 



/Jl {iti)"'{t3-ti){ti-t2) 



(42) 



A more detailed description of the higher-order Bessoid 
integrals as well as the derivation of the recurrence re- 
lation and the amplitude equations can be found in ap- 
pendix^ 
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FIG. 3: Refraction of a ray — propagating from Q to P — 
by a sphere with radius a and refractive index n. The picture 
is drawn in the meridional plane and all indicated angles are 
positive. 



IV. THE SPHERE 
A. Geometrical optics solution 

Consider a plane wave falling on a transparent sphere 
in vacuum. Figure 13 illustrates the refraction of a sin- 
gle ray in the meridional plane, containing the point of 
observation P and the axis. Within the frame of geomet- 
rical optics the cuspoid is formed behind the sphere in 
analogy to figure |21 

Let a be the sphere radius and rt > 1 its refractive 
index. In contrast to the previous sections, we choose 
the origin of the axially symmetric cylindrical coordinate 
system (p, z) differently now, namely as the center M of 
the sphere. The incident plane wave propagates parallel 
to the z-axis. The geometrical optics focus, formed by 
the paraxial rays, is located at = (0, /) with- -^ 



/ 



2 n - 1 



(43) 



A ray passes the point Q, is first refracted at Qi, a sec- 
ond time at Q2 and propagates to P. The incident and 
transmitted angle, 9i and 9t, are related by Snell's law, 
sin^j = nsin^t. Writing the position of P = {p, z) in 
polar coordinates, p = I sinO and z — I cos 9, one can find 
the following expression, determining the three rays that 
arrive at P: 



I sinf 



2 9.1 -2 9t) =asm( 



(44) 



where one has to substitute 9t = arcsin[(sin0i)/n]. This 
is a transcendental cubic-like equation which has three 
roots, either all real or one real and two complex conju- 
gate. (For n > \/2 this is true for z > a; if n < \/2, 
the 3-ray region does not start until some distance be- 
hind the sphere.) We denote them as 6'ij (j — 1, 2, 3) and 
choose their order consistently with the previous nota- 
tions. Therefore, 9i^\ is always real and negative, whereas 
9i^i and ^^^3 are either real and positive (lit region) with 
^i,2 < ^i,3 or complex conjugate (geometrical shadow). 

When the 9i^j are known, we find the 9t.2 from Snell's 
law and the a, and from 



(i^29. 



(45) 



Omitting the index j, the three ray coordinates can be 
written as 



QiP 



I cos( 



a cos OL 



cos /3 



(46) 



The eikonal is the optical path accumulated from Q to 
P (on the dashed vertical line in figure |31 all rays are still 
in phase): 



i> = QQ\ + n Q1Q2 + Q2P - a 
= a{2n cos 9t — cos 9i) + s . 



(47) 



The sphere radius a was subtracted from the path con- 
tributions to make the eikonal zero in the center M, if 
there were no sphere. 

Next we calculate the geometrical optics amplitudes 
by determining the meridional and sagittal radii of cur- 
vature, Rm and Rg, and their changes due to refrac- 
tion. Formulas for the refraction on an arbitrary surface 
with arbitrary orientation of the main radii exist in the 
literaturei^^*^ A simple derivation for the sphere can be 
found in appendix ^ It yields the dependence of the 
actual radii of curvature Rm and Rg (right after the re- 
fraction) on the initial radii Rmo and Rso (just before the 
refraction) : 



Rm 

Rs 



n a RmQ cos^ 9t 



a cos^ 9i + R„io (cos 9i — n cos 9t) ' 

n a Rso 

a + Rso (cos 9i — n cos 9t) 



(48) 
(49) 



For a plane wave, i?mOi RsO ~* 00, the radii of curvature 
in the points Qi (inside the sphere) and Q2 (outside the 
sphere) have the compact form 



R„ 



R-n 



R. 



sm Ui cos Ot 

sin(0, -9t)' 

cos 6i ( cos 9i sin 9t 



2 \sm{9^-9t) 
sin 9i 



R 



S,Q2 



sin(0, - 9t) ' 
sm{2 9t~9^) 
sm{2 9,- 2 9t) 



(50) 

1), (51) 
(52) 
(53) 



The overall geometrical generalized divergence after both 
refractions reads (index j omitted) 



1 



m,Qi 



VJ ^/{Rm.Q^ + d) {Rs^Q, + d) 



V(Pm,Q2 + s) {Rs,Q2 + S) 



(54) 



where d = 2a cos 9t is the distance of propagation within 
the sphere. Note that ray 1 has a negative angle 9i. 
Besides, a double caustic phase shift should be added 
(manually) to the phase of this ray (minus sign) as in 
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FIG. 4: (a) A ray propagates from Q to P in the merid- 
ional plane (plane of incidence), (b) Decomposition of the 
initial electric field vector with length Eo into its n- and a- 
component parallel and perpendicular to the meridional plane 



(|18|l . The caustic shifts of the rays 2 and 3 are taken into 
account automatically if the branch cut for the square 
roots in H54|l is along the negative real axis from —oo 
to and the branch with = +i is used. In this 

procedure it is not allowed to multiply the radicands and 
write them under one common square root. Also the 
case of complex rays 2 and 3 is covered correctly by this 
convention. 

Finally, the geometrical optics solution for the sphere 
is given by (|2()(l . where the eikonal ^p and divergence J 
are given by (|47|l and H54|l . The equation determining 
the three rays is H44() . In the geometrical shadow the 
sum ()20|l becomes only the term with j — I. 

To incorporate Fresnel transmission coefficients, we as- 
sume that the incident light is linearly polarized in x- 
direction, i.e., the incident electric field vector is 



(55) 



with the unit vector in x-direction and Eq = Uq. Since 
axial symmetry is broken, we introduce the polar angle tp 
which is measured from x to y. The point of observation 
P = (p, if, z) will be reached by three rays (two may 
be complex) and their angles Oij are still determined 
by 144|l . for all three rays lie in the meridional plane, 
containing P and the z-axis (figure The initial tt- 
and tr-polarized components depend on ip (figure 01d): 



Eo^n = Eq cos (fi , 
Eo^a = Eq sin if . 

We define the overall transmission coefficients 



12, <T ^21.(7 



= l-r 



12, CT 



(56) 
(57) 



(58) 
(59) 



Here the ti2 ("^12) are the standard Fresnel transmission 
(refiection) coefRcienta^ from the medium 1, i.e., vac- 
uum, into the medium 2, i.e., the sphere. 

The ray field behind the sphere is found by the pro- 
jection onto the original Cartesian system {x,y,z). We 



write the components of the transmission vector T = 
{Tx, Ty, Tz) and show the ray index j = 1, 2, 3 explicitly. 
The (/^-dependence is indicated with the superscript (m) : 



T, 



(0) 



^(2) „ T^(0) 

- T^' ' cos 2 , 



sin2(^. 



T, 



(1) 



Tjrj cos/3,- + Ta 



T^j sin/3,- , 



Tz,j = rj^^ cos(^. 



(2) _ Trr,j cos/3.,- — Tg^j 



(60) 

Hence, the geometrical optics solution for the electric 
field E = {Ex,Ey,Ez) — including the eikonal ^ EZI) 
and divergence J H54II — reads 



Ex 



E,, 



E 



(0) 



E^^ cos 2 (/3 , 



E^^ sin2(p, 
E] ' cos v? , 



(m) 



Eq^ 



(61) 



B. The Bessoid matching solution 




Matching each term = Y^j=\ by the Ansatz 
H21|l in its higher-order formulation 14()|l with the cor- 
responding integral we obtain the vectorial electric 
field E = {Ex, Ey, Ez) in the form 



E = 



(62) 

FigureElillustrates the intensity, i.e., the absolute square 
of the electric field \E\'^ = EE*, for ip = (a;,z-plane) 
and if — t:/2 (y,2;-plane). 

The magnetic field H can be calculated similarly (in- 
cident magnetic field Ho = H^By, Hq — Eq) and the 
(normalized) Poynting vector is given by S = Re(ExH*). 



C. On the axis 

On the axis the electric field is given by its x- 
component only (direction of polarization) due to aver- 
aging over ip in H62|l. For z < f (inside the cusp) it 
is given by the analytical expression H37|) . After several 
simplifications^^ it can be written as 



E = Eq 



TiDi 



i/-l 
Z 



To 



1 - z/f 



(63) 



where the transmission factors Tj = tJ"' are given in 
(|60|l and for dielectric spheres have the form: 

n (1 + 3 cos/3i) cos6'i^i cosSj^i 



To 



{n cos ( 



cos 6t,iY 



An 



(64) 
(65) 
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FIG. 5: Normalized intensity \E/Eof in the normalized x,z- 
plane (top) and in the i/,z-plane (bottom). Contour shadings 
go from white (zero) to black (« 700). Parameters: refractive 
index n = 1.5, dimensionless wavenumber ka = 100. The 
initial electric field vector is Eq = EqBx. The sphere with 
radius a is situated in the origin. The focus of geometrical 
optics is located at 2: = / = 1.5 a, whereas the diffraction 
focus (the point of maximum intensity) is significantly shifted 
towards the sphere: fd ~ 1.25 a. In dimensional units, for a 
wavelength of A = 0.248 the sphere radius is a ~ 4 ^m. 



The phases in the coordinate Z = —2 ^/fi — (p2 are 



fi = tp3 = kai2n cos 9t,i — cos 9i^i + '^^ ) , (66) 



sin /3i 



ip2 — '2 k a (n — 1) + k z . 



(67) 



and Di = \J^\ — ^zj \/J\ is the first ray's compensated 
sagittal divergence: 



\/{Rm.,Qi)l {Rs,Qi)l 



\/[{Rn.,Q,)l+di\[{Rs,Q,)l+di] 
\/{Rm,Q2)l {Rs,Q2)l 



V(^nM3Jl+Sl 



2 fc sin /3 



-2 V2fca cos(/3i/2) 



cot 9i^i cos 6*4^1 sin(/3i/2) 
1 + l/n2 - 3sin2(/3i/2)/ sin^ 6*,,! 



, (68) 




|£(/d)P/£o 
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FIG. 6: Left: Diffraction focus in units of the sphere radius as 
a function of n and k a (contour lines from top to bottom go 
from 1.1 to 3.0 in steps of 0.1). Right: Intensity enhancement 
at fd (contour fines from bottom left to top right are 20, 50, 
100, 200, 500, 1000 and 2000) 



which manifestly has no singularity until the geometri- 
cal focus where iRm,Q2)i + Si ^ 0. The minus sign 
comes from the manually inserted phase shift of the first 
ray. (All aforementioned quantities should be expressed 
in terms of the angles 9i^i and 9i^2 = as described in 
detail in section llV DI below.) The structure of the first 
two lines in H68(l is general and is valid for arbitrary axi- 
ally symmetric systems. Di is always finite on the axis, 
since both the sagittal radius of curvature and the phase 
difference ipi — are proportional to the distance p. 

Equation H63II is valid even near the focus, since the 
diverging terms Di/Z and (1 — z/f)^^ almost cancel. 
For z ^ f, however, the divergence of Di itself becomes 
important, as the non-paraxial ray 1 becomes axial. 

In figure El we show the position and the value of the 
maximum of \E\^ as a function of the refractive index 
and the dimensionless product ka, calculated from H63|l . 
The z-coordinate of this global maximum is denoted with 
fd {diffraction focus) and the intensity there is \E{fd)\ ■ 
Contrary to the square dependence for the case of an 
ideal lens^ even for macroscopic spheres the maximum 
intensity turns out to be about proportional to k a, in 
agreement with the general theory. 

The main contribution in (|63|l stems from the Bessoid 
integral, that is from the term oc Ti Di I. Thus, the po- 
sition of the maximum can be estimated from condition 
(|35|l . i.e., — </?2 ~ 3 7r/4. If the phase difference ipi — ipi 
from H65|l and (|67|l is expressed as a function of 9i x^ Tay- 
lor expanded and equated to 3 7r/4, then we get in the 
lowest non-trivial order of the inverse product k a: 



3 vr n (3 — ri) — 1 



2 n- 1 



4 fc a n{n — I) 



(69) 



Hence, in the limit of small wavelengths or large spheres 
the relative difference between the diffraction and the 
geometrical focus decreases proportionally to the inverse 
square root of ka. The factor 3 7r/4 « 2.356 can be re- 
placed by the more exact Bessoid value 2.327 from H35|l . 
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Expression (|69() approximates the position of the maxi- 
mum within an error of < 5 % for k a > 100 and values of 
the refractive index in the range 1.4 < n < 1.6. The tran- 
scendental phase difference condition (|35l) . which holds 
for large angles, naturally has a wider range of appli- 
cability. With very good accuracy the diffraction focus 
also provides the maximum for the absolute square of 
the magnetic field, \H\^ = HH*, as well as for the z- 
component of the Poynting vector S. Note that on the 
axis H — H By and S = Se^- 

D. A protocol for the electromagnetic field 
calculation behind the sphere 

For convenience we summarize the sequence of steps 
that should be used for the calculation of the field behind 
a sphere irradiated with linearly polarized light on the 
basis of the formulas developed above: 

1. Finding the rays: The origin of the coordinate 
system is in the center of the sphere. Choose a 
point P = (/o, if, z) behind the sphere and numer- 
ically calculate the 3 rays arriving at P. These 
rays are characterized by the 3 incident angles Oij 
(j = 1,2,3), found numerically from 144|l . and 
numbered according to figure El AH other angles 
follow from Snell's law and from (|45|l . respectively. 
Outside the cuspoid the rays 2 and 3 are complex. 

2. The geometrical optics solution: Compute the ge- 
ometrical optics solution for the electric field (|61|l . 
which is the sum of contributions from three rays. 
The eikonals ipj are calculated from (|47|) with (|46|l . 
The geometrical optics amplitudes are given by the 
Fresnel transmission components tJ™-* (I6()|I and the 
generahzed divergence factors l/-\/J^ itK^ . which 
follow from the radii of curvature H50|) - (|53|l and the 
distances Sj from the sphere to P (^BJ. The con- 
ventions for the complex roots shall correctly add 
up all individual caustic phase shifts. Ignore the 
fact that the geometrical field diverges near caustic 
regions. 

3. Bessoid matching: Starting from the eikonals ipj, 
determine the Bessoid coordinates, i.e., first Z, and 
then R and x H28|l . Next, compute and correctly or- 
der the points of stationary phase tj (126(1 , most con- 
veniently using trigonometric formulas.^— With the 
generalized divergences (|54|) and the Hessians (|25f) 
the Bessoid amplitudes Am, AmR and A^z can 
be computed from (|42|l for all orders m = 0,1,2. 
The electric field component iJ^™) associated with 
the m-th order results from the Ansatz H40|l . The 
Bessoid-matched field E is finally given by (|62|l . In 
the case of a scalar plane wave (or exactly on the 
axis) only the order m = contributes. Proceed 
accordingly for the magnetic field H, employing dif- 
ferent transmission coefficients in step 2. 




1.2 1.4 1.6 1.2 1.4 1.6 

FIG. 7: \E/Eof on the axis. Dashed lines represent the Mie 
theory, solid lines are the results of Bessoid matching. The 
parameters are as in figure |S] and the cases (a), (b), (c) and 
(d) correspond to q = 300, 100, 30 and 10, respectively. In di- 
mensional units, for A = 0.248 fim this corresponds to sphere 
radii of a ~ 12 /im, 4/im, 1.2 /im and 0.4 /im. 

4. The Bessoid integral: The final solution H62|) con- 
tains the Bessoid integral and higher-order Bessoid 
integrals as well as their partial derivatives. The 
Bessoid integral I{R, Z) can be efficiently com- 
puted numerically via the differential equation l|12|l . 
The higher-order Bessoid integrals follow from a re- 
cursive relation H39|l . 

5. Remarks: The speed limiting bottleneck of this 
procedure is the finding of the rays 9ij in step 
1. The numerical evaluation of the Bessoid inte- 
gral is very efficient and in all other steps analyti- 
cal expressions are applied. When approaching the 
caustic (axis, cuspoid) , individual quantities — the 
geometrical optics amplitudes — diverge but their 
combinations remain finite. In our calculations we 
observed perfect numerical stability up to distances 
from the caustic of the order of 10"^ times the 
sphere radius, which is more than enough for any 
practical purposes. 



E. Comparison with the theory of Mie 

We presented a general way to match geometrical op- 
tics solutions with the Bessoid integrals. It can be ap- 
plied to any axially symmetric system with the cuspoid 
topology of spherical aberration. 

For the sphere we can compare our approximate results 
with the theory of Mie>i5. A main quantity characterizing 
the sphere is the dimensionless Mie parameter q = ka. 
FigureEJcompares the intensity on the axis obtained from 
the Mie theory with the Bessoid approximation. The 
parameters are as in figure [S] and the Mie parameter is 
q = 300, 100, 30 and 10. 
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FIG. 8: Normalized \Ef, \H\^ and Sz in the normalized x,y- 
plane for z = 1.02 a calculated with Bessoid matching (left) 
and with the Mie theory (right). The parameters are the same 
as in figure |S] 



We see very good agreement down to 5 « 30 (a/ A « 
4.8). For q = 10 (a/A « 1.6) the asymptotic behavior far 
from the sphere is still correct. However, for small q the 
characteristic scale a is no longer large compared to the 
wavelength A and geometrical optics becomes invalid. 

Next, we compare the off-axis electric and magnetic 
field as well as the z-component of the Poynting vector, 
Sz (figure ISJ. Right behind the sphere (z — a) the agree- 
ment is not perfect (see figure [Tj), though all qualitative 
features are preserved. Sections at z = 1.02 a already 
show good agreement (figure |SJ) and for z > 1.05 a the 
pictures become visually almost indistinguishable. Dis- 
cussing the quality of these results, one has to differenti- 
ate between the accuracy of the method and the influence 
of those factors which can be taken into account, but were 
not included into the current consideration. 

The accuracy of the Bessoid matching procedure it- 
self was studied separately for the case of a spherically 
aberrated wave incident onto an aperture. The deviation 
— defined as the maximal relative error of the intensity 
(|£'Bossoidp-|£^oxactn/|£^cxactP — bctwccn the Bcssoid 
matched geometrical optics solution and the correspond- 
ing (exact) Rayleigh-Sommerfeld diffraction integraliSi 
decreases as the aperture increases. If the aperture is 
large enough, the deviation is below 10"'^ for spherical 
abberation strengths and wavevectors approximately cor- 
responding to the focusing by spheres studied in figures 



□ andil 

For the sphere, for all investigated Mie parameters 
30 < (7 < 300, the deviations from the exact Mie solution 
are about ±5 % in those regions where the intensity is not 
very small (including the caustic axis and cuspoid). A 
detailed analysis indicates that this deviation originates 
from several factors: 

(i) Influence of the finite size of the sphere. There exist 
diffractive contributions from creeping rays^ propagat- 
ing along the sphere surface. They can in principle be 
accounted for at the expense of the simplicity of the pro- 
cedure, essentially by considering the interference of the 
Bessoid field with these additional rays. This results in 
oscillations which can actually be seen in the Mie curve 
on the right side of the Bessoid tail in figures [7^ and b. 

(ii) Rays entering the sphere undergo multiple inte- 
rior reflections. Some of them satisfy resonance con- 
ditions, accumulate significant energy inside the sphere 
and refract outside. This produces an additional field 
behind the sphere, in particular the intensity becomes 
non-zero in the regions of geometrical shadow for the 
directly transmitted rays used in the Bessoid matching. 
Here again, one can (in principle) study such multiply re- 
flected rays separately and add them to the Bessoid field, 
but their contribution to the focusing properties of the 
microspheres is of secondary importance. 

(iii) Finally, very close to the sphere surface at dis- 
tances of the order of a fraction of A, there exist evanes- 
cent contributions, which are taken into account in the 
Mie theory, but are obviously absent in the Bessoid 
matching procedure. 

Thus, the quality of the Bessoid matching in the most 
interesting regions near caustic surfaces is quite satis- 
factory. It rectifies the divergences of geometrical op- 
tics, which is asymptotically correct for large q in non- 
singular regions of space. Clearly, the procedure has to 
be extended whenever contributions from rays other than 
those 3 used for matching become significant. 

The field distribution behind the sphere has a rich fine 
structure (figure|Sl) which our geometrical approach helps 
to clarify. It is known that the ring-type field enhance- 
ment corresponds to the cuspoid caustic, having the ap- 
proximate radial distancaU 

(4-n2)3/2 



3^/3? 



(70) 



Our approach explains the double-peak structure of 
along the direction of polarization. It is related to the ax- 
ial field component and can be understood in terms 
of geometrical optics. On the axis, the Ez components 
from the rays 1 and 3 point into opposite directions and 
cancel, having an effective phase difference of tt. Off the 
axis, ray 1 underwent a caustic phase shift of — 7r/2 when 
crossing the axis, which makes the condition for construc- 
tive interference: ipi — 953 = 37r/2. Then, according to 
(|34|l . the peak occurs at the radial distance 

_ - f3 ^ 3 A 
^P'"2fcsin/3 8 sin/3' ^ ^ 



12 



More details on the derivation are given elsewhere^ to- 
gether with the refined coefficient 0.293 (instead of 3/8) 
obtained from the Bessoid asymptotic ljE8(l . 

Double-peak structures have been observed in 
nano-patterning experimentsS*^2iiS and were semi- 
quantitatively explained on the basis of the Mie solutiontS, 
In an actual experiment it may depend on the laser pulse 
parameters and the properties of the patterned material, 
whether the Poynting vector or the electric field is re- 
sponsible for the patterning process. For small spheres 
this double peak effect can be understood using the near 
field pattern for a scattering dipole)S*iS The present ex- 
planation (for sphere diameters of a few wavelengths and 
larger) results in the same orientation of the maxima and 
thus these two limiting cases cover almost all range of 
sphere sizes. Similar polarization dependence of the field 
distribution in focal regions can be used to improve the 
resolution4i 



V. CONCLUSIONS 

We described theoretically arbitrary axially symmetric 
aberrated focusing and studied light focusing by micro- 
spheres as an example. Following the method of uniform 
caustic asymptoticSfi^i we introduced a canonical inte- 
gral describing the wave field for the given cuspoid ray 
topology. This Bessoid integral appears naturally in the 
paraxial approximation. In some regions (off the caustic 
or exactly on the axis) it reduces to simple analytical ex- 
pressions. In other regions we efficiently computed this 
highly oscillatory integral via a single ordinary differen- 
tial equation. 

For arbitrary axially symmetric focusing, coordinate 
and amplitude transformations match the Bessoid wave 
field and the solution of geometrical optics. The caustic 
divergences of the latter are removed thereby. For vecto- 
rial problems with angularly dependent field components, 
higher-order Bessoid integrals are used for the matching 
procedure. The formulas significantly simplify on and 
near the axis. An approximate universal condition for 
the diffraction focus can be given in terms of phase dif- 
ferences. Here, the concept of caustic phase shifts is of 
main importance. 

The central part of the Bessoid integral is essentially 
a Bessel beamSS, with a variable cross section due to the 
variable angle of the non-paraxial rays. Its local diameter 
is always smaller than in the focus of an ideal lens with 
the same numerical aperture. Besides, the largest possi- 
ble apertures can be physically realized, which is hardly 
possible with lenses. All this is achieved at the expense 
of longitudinal confinement. 

As an example the focusing of a linearly polarized 
plane wave by a transparent sphere is studied in de- 
tail. We calculate the geometrical optics eikonals and 
divergences, incorporate Fresnel transmission coefficients 
and perform Bessoid matching. Using the general theory, 
simple expressions for the light field on the axis and for 



the diffraction focus are derived. The two strong maxima 
in the intensity observed immediately behind the sphere 
can be explained as well. 

Finally, the results of the Bessoid matching procedure 
are compared with the Mie theory. The agreement is 
good for Mie parameters ka > 30. Near the sphere the 
correspondence is worse due to unaccounted evanescent 
contributions. 

The developed formulas can be directly applied in 
other areas of physics where non-paraxial axially sym- 
metric focusing is of importance, e.g., acoustics, semi- 
classical quantum mechanics,'*^ fiat superlenses based on 
left-handed materials)^ radio wave propagation, scatter- 
ing theory^Si chiral conical diffraction^ etc. 

Concluding, let us briefly enumerate several possibili- 
ties to extend and reflne the developed formalism. Weak 
absorption can be incorporated easily, for it just changes 
the amplitudes along the rays and the transmission coef- 
ficients. Strong absorption additionally modifies Snell's 
law of refraction, still preserving the axial symmetry. 
One can consider incoming radially or azimuthally po- 
larized beams, which are known to produce better reso- 
lution than linear polarization^! The diffraction of light 
from regions beyond the sphere radius can be incorpo- 
rated by considering the interference of the Bessoid field 
with creeping raySf22i For other geometries, in particular 
finite apertures with sharp boundaries, edge rays or the 
Rubinowitz representation^ or an approach based on 
catastrophe theory^ have to be used. Such corrections 
become relevant, for example, for the ray structure and 
the field distribution immediately behind spheres with a 
refractive index n < \/2. Finally, one can calculate the 
interference of the diffracted light with the original in- 
cident wave or the interference of the light refracted by 
several spheres or arrays of spheres. The latter yields 
interesting secondary patternaiSi related to the so called 
Talbot effectii 
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APPENDIX A: A NEAR AXIS APPROXIMATION 
FOR THE BESSOID INTEGRAL 

We make the substitution w = p\ in ||2J): 

1 Z — + — 
I{R,Z)^- Jo(i?V^)e V 2 ^J^yj_ (Al) 

2 Jo 

Near the axis (smaU R) the Bessel function is slowly vary- 
ing compared with the exponent. The integral will have 
significant contribution only from the region in which the 
exponent's phase is stationary, i.e., regions near w = —Z. 
We consider the most interesting caustic part of the axis 
for which Z < 0. In a lowest order approximation the 
Bessel function is considered as constant near the sta- 
tionary point —Z and can be pulled out of the integral. 
The phase can be written as a complete quadratic form. 
With the full square of v = {w + Z)/2: 



I{R,Z) w Jo{R^/^) 



4 



dw . 



Z/2 



(A2) 



The remaining integral can be expressed in terms of the 
complementary error function^'^ erfc (of complex argu- 
ment) and hence we arrive at (|ll|l . 

For Z > the point w ~ Q should be taken as a sta- 
tionary edge point of the integration. And the near axis 
approximation Ijlll) remains valid as long as the Bessel 
function is set to Jo(0) = 1. 



APPENDIX B: AN ORDINARY DIFFERENTIAL 
EQUATION FOR THE BESSOID INTEGRAL 

We derive the paraxial Hclmholtz equation 



Irr + ^Ir + 2iIz 



(Bl) 



as well as the following ordinary differential equation for 
the Bessoid integral: 

lBRR + ^lRR-(^-^ + Z^lR + iRI = 0. (B2) 

Indices denote partial derivatives. Both equations can be 
rewritten in the compact form 



L + 2iIz - 0, 
Lr - ZIr + iRI = 0, 

where L is the radial Laplacian 
i = Irr 



(B3) 
(B4) 



(B5) 



We begin with the proof of (|Bip and state that we 
may differentiate under the integral sign, since the par- 
tial derivatives of the integrand exist and are continuous 



functions. Starting from the Bessoid integral in the polar 
representation its integrand can be written as 



G = piMRpi)E, 
with the abbreviation 

/ 2 4 



E = c 

The (multiple) partial derivatives are 
GR^--plJiiRpi)E, 



G 



RR 



[MRpi)-J2iRpi)]E, 



pI 



Gz = -i^JoiRpi)E. 



(B6) 



(B7) 



(B8) 
(B9) 



(BIO) 



Here we used the derivative formula for Bessel functionsS^ 

(Bll) 



d /.N _ Jm-l(f) — Jm+l{t) 



with m = to obtain 1B8|I and m = 1 for (|B9|I . Note 
that J-i[t) = ~Ji{t). Applying the recurrence relation 
for Bessel functionsS^ 



2 m 

Jm+l{t) = ^Jm-l{t) H — Jm{t) : 



(B12) 



one can eliminate J2 from (|B9|I . And then it is enough 
to notice and verify that 

GRR + ^GR + 2iGz^0. (B13) 

This proves equation IjBip . 

For the proof of ljB2|l . we need to note that its left 
hand side can be expressed as the integral of a partial 
derivative 



H 



d 

dpi 



\pi Ji{Rpi) e 



^ 9 + /I 



dpi 



(B14) 

With the help of (|BTT|) and (jBT2|) both the left hand side 
of (IB2II and H become 



/•oo 

/ [iRpiJQ{Rpi) + pl{Z + pl),h{Rpi)]Edpi. (B15) 



Thus, in order to prove (|B2I) . it is enough to show that 
H = Q. This follows from the Newton-Leibniz formula 
applied to the (definite) integral ljB14p : 



H 



ipi Ji{Rpi) e 



(B16) 



The lower bound at vanishes for obvious reasons. For 
the upper bound at 00 one assumes an infinitely small 
imaginary part in front of the fourth order term in the 
exponent: p\ {l—ie)p\ with e > 0. This completes 
the proof of the differential equations (jBip and ljB2p for 
the Bessoid integral. 
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and thus: 



np 



FIG. 9: Near the axis, the phases of the rays 1 and 3 differ 
from the phase of the non-paraxial rays np (crossing the axis 
at the same z) by ±k p sin /3, whereas the phases of ray 2 
and the (par) axial ray are the same in first order. The non- 
paraxial ray and ray 3 cross the axis at an angle /3 > 0. 



APPENDIX C: THE NEAR AXIS BESSOID 
COORDINATES 

Near the axis the phases of the rays can be Taylor 
expanded. From figure El one infers that up to the first 
order in p tlie phases can be written as 



ifii ~ Vnp + fcp sin/3, 
V's ~ y^np ~kp sin f3. 



(CI) 



Here ip^p and ipp denote the phases of the non-paraxial 
rays and the (par)axial ray (with p = 0) and /3 > is the 
angle of ray 3 with the axis. 

We insert these phases into the exact expressions for 
R and Z in (|28|l . Taylor expand the result with respect 
to p and resubstitute ipnp « (vi + '/'3)/2, (Pp ~ V2 and 
kp sin 13 {(fii — (p3)/2 from HC1|I . This yields (|22J) and 



APPENDIX D: THE ON AXIS FIELD 

Here we derive a simple on axis expression for the 
Bessoid-matched field U (|21|1 . namely equation H37I) . 

On the axis and inside the cusp p = (i? = 0) and 
2r < (.^ < 0). The stationary points, given by (|26|l . are 



ti = -V^, <2=0, t3 = -ti. 

Then, the amplitude A in (|31|l simplifies to 



A = C/o ■ 



/J2 



(Dl) 



(D2) 



because on the axis the ratios •\/ -ffi,3/-\/t/i.3 are both 
finite and the corresponding other terms disappear upon 
multiplication with t2 = 0. Due to the restriction to the 
lit region {Z < 0), all rays are real and equation H24|) 
holds. By virtue of Q detH2 = Z'^, and due to l(TUjl 
signH2 = 2, one finds 



iZ . 



(D3) 



A = 



Ua.2Z 



(D4) 



This approximation for the amplitude A is valid up to 
the focus [Z ~ 0). As ray 2 converges like the inverse 
distance from the focus, \/J2 is proportional to Z . 

The amplitude Aji in H31II vanishes due to = —ti 
and the fact that 



(D5) 



which means that rays 1 and 3 have equal amplitudes 
and that the caustic phase shifts are in accordance with 
the signature of the Hessian. Consequently, 



Ar = 0. 

With HD1|) and (jPSp the amplitude Az reads 



Az = -\U, 



^04 



0.2 



H2\ 



J2J 



(D6) 



(D7) 



The first term is non-trivial. Both \/ Hi and ^fJi are 
zero on the axis, but their ratio is finite and well defined. 
Indeed, the Taylor expansion of Cardan's solution ti in 
its trigonometric representation^^ yields in the first order 
in R: 



(D8) 



Therefore, again in first order in R: det Hi 
Due to sign Hi = —2 we obtain 



'Hi = i V2-R 



and with (ID 3 



A^^2 



Uo,iV2R^ 



ZVJi 



--^-2i^ 



(D9) 



(DIO) 



This approximation for Az holds for small values of R. 
It is finite, since \/Ji approaches zero as for i? ^ 0. 

For the final representation of the field U, we can sub- 
stitute the near axis expression for R-\J —Z 134|l into Az ■ 
On the axis the phase coordinate becomes X — ^2i which 
results from substituting </?i — (^3 and Z = —2 \J^\ — ^2 
into the corresponding expression in (|28|l . This leads to 



Al^ 
^/T2 



■ Ar Ir 



{iZI-2Iz 



Aziz 



(Dll) 



2 Uo,i V'^kp sin f3 
Zy/J[ 



Using the linear relationship between the Bessoid integral 
and its Z-derivative ((T3|l . 



iZI-2Iz = 1, 
we end up with equation H37|) . 



(D12) 
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APPENDIX E: HIGHER-ORDER BESSOID 
INTEGRALS 

Higher-order Bessoid integrals (|38|l appear naturally, 
if one expands an arbitrary initial field amplitude on the 
aperture in a Fourier series: 

oo 

Uo{pi,fi) = ^ [am(Pi)cos(m^i) + 6„(pi)sin(mi^i)] . 

m— 

(El) 

The form of the coefficients am and bm can be seen 
from a two-dimensional Taylor expansion in Cartesian 
coordinates around the point (0,0), rewritten into polar 
coordinates: 



Uo{pi,fl) = 'Y Cran pT cos"" "lfiiSm"lfii 

ni—O n—0 

with 



1 fm\ 9"t/o(4,y;) 



(E2) 



(E3) 



x[=0, y[=0 



Thus, is the lowest possible power of pi which can 
be found in the term with expiirmpi). An additional pi 
comes from the transformation from Cartesian to polar 
coordinates. 

If we define the functions 

/:„,EE/„e'™'^, (E4) 

we find that they satisfy the paraxial Helmholtz equation 

(E5) 



2 i Im,Z + ImMR + Im,B. + Im,ipip — , 



where = Irn- 

Due to HB11|I and (|ljl2() . one can write the identity 

PT+' Jrn+l{Rpi) = -^[Pr' JmiRpi)] 

+ ^pT+'j,n{Rpi). (E6) 

Hence, the recursive relation for the Bessoid integrals 
follows: 



(E7) 



i-e., Ii = - Io,R = -Ir, h = Irr-Ir/R, etc. Using ^7} 
and ljB12|l as well as expression (|ll|l for /, one obtains 

Im{R, Z) K, ^ Jm{RV — Z) 



2 

■ Z^-^ / Z ■ 2L 

X e' 4 erfc — e' 4 



(E8) 



which is the analytic near axis expression for the higher- 
order Bessoid integrals. 




FIG. 10: Meridional cross section for the determination of the 
meridional radius of curvature, Rm- M is the center of the 
sphere 



While the coordinates and phases (i?, Z, x) remain 
unchanged, the derivation of the higher-order amplitudes 
(v4„i, AmR, Amz) rcquircs some insight for m > 2. Let 
us briefly consider the case to = 2. For the matching 
procedure we need the asymptotic behavior of I2 = Irr — 
Ir/R far from the caustic regions where i? ^ 1 and 
where it is dominated by the term Irr. Note that for 
the matching procedure we need exactly this asymptotic 
representation and in the non-caustic regions only. Thus 
— although we need the second-order Bessoid integral I2 
on and near the axis, where it vanishes — we shall use 
its asymptotic stationary phase expressions far from the 
axis for the derivation of the amplitudes. In this region 
it is equivalent to the asymptotic of Irr. 

In fact, we may generalize this statement to arbitrary 
order. Due to l|E7|l the leading term in the stationary 
phase calculation is always 



-—] I- 

dR 



(E9) 



Therefore, the equations for the amplitudes H30|) become 



U, 



(m) 



(ii,r 



A 



m AjjiR 



1+2/1 

2 "-j ^mZ 



(ElO) 



which can be seen from the Bessoid integral's Cartesian 
representation with the phase Q : 

The equations (|E10() have the same form as (|30|l except 
an additional factor (iij)™ on the left hand side, proving 

dm. 



APPENDIX F: WAVEFRONT RADII OF 
CURVATURE FOR THE REFRACTION ON A 
SPHERE 

Let us consider a point source G and start with the 
derivation of the meridional radius of curvature (figure 

cnii. 
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We express 6 from (jF2|) , substitute it into ljFl|l and finally 
obtain the meridional radius of curvature (|48|l 



FIG. 11: Meridional cross section for the determination of the 
sagittal radius of curvature, Ra 



The initial radius of curvature is Rmo = GE, the one 
after refraction is Rm = ED. The infinitesimally neigh- 
bored beam (7^1) which is refracted in E' (the angles 
of incidence and transmission in E are 9i and 9t, in E' 
they be denoted 0- and 9'^) also propagates to D. The 
normals onto Rmo through E and onto Rm through E' 
are g and d, respectively. In the necessary order the 
length of the arc EE' can be approximated by the dis- 
tance e w EE'. As all angles are small: g = j RmO, 
d = 5 Rm, and e = ea. On the other hand, we find from 
the infinitesimal triangles: g — e cos 9i, d = e cos 9t. This 
leads to 



Rn 



7 Rmo cos 9t 
S cos 9i 



(Fl) 



where we have introduced a minus sign because the wave 
is converging after the refraction. The remaining problem 
is the angle S in the denominator. To find 5 we write the 
relations between angles and primed angles: 



%=9,+-/ + e, 



6 + e. 



(F2) 



With Snell's law 



. sin6'i . sin(6', +7 + e) 
arcsm arcsm , (1 3) 



n n 
a first order Taylor expansion in (7 + e) yields 

, (7 + e) cos 9, 
n cos 9t 



(F4) 



Rtn. — 



n a RmO cos^ 9t 



a cos^ 9i + RmO (cos 9i — n cos 9t) 



(F5) 



For the sagittal radius of curvature we consider figure 

A ray which emerged from G is refracted in E (incident 
angle 9i and transmitted angle 9t). The distance from E 
to the intersection H of the ray with the line passing 
through G and the sphere center M is the sagittal radius 
of curvature, as a neighbored ray, emerging from G and 
hitting the sphere not in E but infinitesimally shifted 
perpendicular to the meridional plane, will also propa- 
gate to H due to symmetry around the line GM. We 
have Rso = GE and Rs = EH. The tangent theorem 
states (all angles are in general large now) 



1^ - ^ a - Rso ^ TT - 
tan — - — = — — cot — — - 



a + Rso 



(F6) 



where tt — 9i is just the third angle in the triangle GME. 
Due to v + ^J, = 9i wc find 



a — RsQ TT — I 

Li — arctan cot 

2 \a + Rso 2 



In the triangle MHE the sine theorem reads 

Rs sin(7r — /i) 
a sin 7/ ' 



(F7) 



(F8) 



where we have again introduced a minus sign due to the 
convergence of the refracted wave. Trigonometric trans- 
formations finally give the sagittal radius of curvature 
631 



Rs 



n a Rso 



Rso (cos 9i — n cos 9t) 



(F9) 
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